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When using Hartree-Fock (HF) trial wave functions in quantum Monte Carlo calculations, one 
faces, in case of HF instabilities, the HF symmetry dilemma in choosing between the symmetry- 
adapted solution of higher HF energy and symmetry-broken solutions of lower HF energies. In this 
work, we have examined the HF symmetry dilemma in hydrogen rings which present singlet instabil- 
ities for sufficiently large rings. We have found that the symmetry-adapted HF wave function gives 
a lower energy both in variational Monte Carlo and in fixed-node diffusion Monte Carlo. This indi- 
cates that the symmetry-adapted wave function has more accurate nodes than the symmetry-broken 
wave functions, and thus suggests that spatial symmetry is an important criterion for selecting good 
trial wave functions. 
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I. INTRODUCTION 

It is well known that the symmetry-adapted solution 
of the nonlinear Hartree-Fock (HF) equations of an elec- 
tronic system is sometimes unstable. An unstable so- 
lution corresponds to a saddle point of the energy as a 
function of the orbital parameters, and breaking of space 
and/or spin symmetries of the wave function then neces- 
sarily leads to one or several lower-energy HF solutions. 
The stability conditions of the HF equations were first 
formulated by Thouless [i[, and the different instabil- 
ities were first categorized by Cizek and Paldus [2|-[3]. 
For closed-shell systems, one may encounter "singlet in- 
stabilities" when only spatial symmetry is broken, and 
"triplet (or nonsinglet) instabilities" when spin symme- 
try is also broken. There is thus a symmetry dilemma [8| 
in choosing between the symmetry-adapted wave func- 
tion of higher HF energy and a symmetry-broken wave 
function of lower HF energy, in particular as a reference 
for a post-Hartree-Fock calculation. 

A clear example is provided by closed-shell hydro- 
gen ring s H4„-|.2 with equal bond lengths Q (see, also, 
Ref. KLOl ). The symmetry-adapted HF solution exhibits 
singlet instabilities for sufficiently large numbers of hy- 
drogen atoms, and one can obtain symmetry-broken HF 
solutions with orbitals localizing on cither the atoms or 
the bonds. However, both M0ller-Plesset perturbation 
theory and linearized coupled cluster doubles theory (also 
called CEPA-0 or DMBPT-oo) give a lower total en- 
ergy when starting from the symmetry-adapted solution 
than when starting from the symmetry-broken solutions, 
which casts doubts on the physical significance of the 
symmetry-broken solutions. Of course, the symmetry 
dilemma would be removed with a full configuration- 
interaction calculation which must give one unique so- 
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lution, independent of the orbitals used. 



Quantum Monte Carlo (QMC) approaches are alterna- 
tives to the traditional quantum chemistry methods |ll| - 
Il3| . The two most commonly used variants are varia- 
tional Monte Carlo (VMC) which simply evaluates the 
energy of a flexible trial wave function by stochastic 
sampling, and fixed-node diffusion Monte Carlo (DMC) 
which improves upon VMC by projecting the trial wave 
function onto the ground state subject to the condition 
that the nodes of the projected wave function are the 
same as those of the trial wave function. For large sys- 
tems, the most common form of the trial wave function 
is a Jastrow factor multiplied by a fixed HF determinant 
(though for small systems one can do much better by re- 
placing the HF determinant by a linear combination of 
optimized Slater determinants). If a system exhibits HF 
instabilities, then QMC also faces the symmetry dilemma 
in choosing between different HF wave functions. Indeed, 
different HF wave functions necessarily lead to different 
energies not only in VMC, but also in DMC since the 
nodes of these HF wave functions are generally differ- 
ent. This symmetry dilemma in DMC is only due to the 
fixed-node approximation, since without this approxima- 
tion DMC would give one unique solution, independent 
of the trial wave function. Of course, for systems that 
are not very large, symmetry breaking could probably 
be avoided in the first place by optimizing the orbitals 
within VMC [11 [11, instead of using fixed HF orbitals. 



In this work, we study the impact of the HF symmetry 
dilemma for QMC in hydrogen rings H4„+2. In Sec. II, 
we review the HF symmetry-breaking problem in these 
systems, and discuss the effect of using a Slater basis 
versus a Gaussian basis. In Sec. HI, we explain the QMC 
methodology and report our VMC and DMC results. Our 
conclusions arc summarized in Sec. IV. 
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FIG. 1: Structure of the one-particle density matrix for the HF symmetry-adapted (SA) and symmetry-broken (SB-AC and 
SB-BC) solutions. 



II. HARTREE-FOCK SYMMETRY BREAKING orbitals ip,{T), 



In previous studies [3, [la] , the electronic structure of 
periodic rings oi An + 2 evenly spaced hydrogen atoms 
(with a fixed distance of t^h-h = 0.74747 A) has been in- 
vestigated. The number of hydrogen atoms is restricted 
to 4n -t- 2 in order to obtain a possible closed-shell single- 
determinant solution with 2n + 1 occupied orbitals. The 
symmetry-adapted HF wave function has a metallic char- 
acter (in the limit of an infinite ring) and can be ex- 
pressed with either delocalized canonical orbitals or lo- 
calized Wannier orbitals (which do not have an exponen- 
tial decay). The canonical orbitals, except for the low- 
est one, are doubly degenerate, and in a minimal basis 
the orbital coefhcients are fixed by the cyclic symme- 
try. Besides the symmetry-adapted (SA) solution, two 
different symmetry-broken HF solutions of lower energy 
can be obtained beyond critical ring sizes, when using 
unit cells of 2 hydrogen atoms. One solution corresponds 
to orbitals localized on hydrogen atoms and is referred 
to as the symmetry-broken atom-centered (SB- AC) solu- 
tion, while the other corresponds to orbitals localized on 
bonds and is referred to as the symmetry-broken bond- 
centered (SB-BC) solution. The SB-BC solution has the 
lowest HF energy. The three solutions can be schemat- 
ically described as • • • H- • • H- • • , • • • H+ • • • H~ • • • , and 
• • • H — H- • • . In each case, the symmetry breaking is ac- 
companied by an opening of an energy gap between occu- 
pied and virtual orbitals, and orbitals decay much more 
rapidly than for the symmetry-adapted soluti on |9| , in 
agreement with the theoretical result of Kohn [13j . 

In order to distinguish the three different wave func- 
tions, one may look at the one-particle density matrix 
P 



PaB = 2 



/ ^ Cq, i Cp i , 



(1) 



v*(?) 



E 



Xa(r), 



(2) 



containing the coefficients Cai of the occupied molecular 



expanded in a minimal set of atom-centered basis func- 
tions, i.e. one single basis function Xq(?) Pcr hydrogen 
atom. As depicted in Figure [1] for the SA solution, we 
see equal elements on the diagonal and the sub-diagonals 
of the density matrix. For the SB- AC solution, an alter- 
nation of element values on the diagonal of the density 
matrix is obtained, but equal elements on the first sub- 
diagonal, and for the SB-BC solution we have equality 
of the diagonal elements and alternation on the first sub- 
diagonal. 

In Ref. y, a minimal Gaussian basis set (five s Gaussian 
functions contracted to one single basis function for each 
hydrogen atom) was used. However, Gaussian basis func- 
tions are not appropriate for all-electron QMC calcula- 
tions. They give large statistical fiuctuations due to their 
incorrect vanishing gradient at the nuclear positions. It is 
thus much preferable to use Slater basis functions which 
correctly have a non-zero gradient on the nuclei and an 
exponential decay at large distance. In this work, we 
use a minimal Slater basis set (one Is Slater function on 
each hydrogen atom) with an exponent of 1.17, which is 
smaller than the optimal exponent of 1.24 for an isolated 
H2 molecule. Spin-restricted HF (and MP2) calculations 
were performed with an experimental code for ring sys- 
tems, employed already for the previous studies |18| . The 
necessary integrals over Slater functions have been calcu- 
lated with the program SMILES 119|. In order to obtain 
the symmetry-broken HF solutions, we start from a set 
of localized Wannier orbitals describing either an ionic 
situation or an explicit bond in the two-atom unit cell, 
and use an iterative configuration interaction procedure 
using singly excited determinants [20, [21| instead of di- 
agonalizing a Fock operator. 

Table [l] reports the HF energy differences between the 



TABLE I: HF energy differences of H4,i+2 rings, per H2 cell in mhartree, between the symmetry-broken solutions (SB- AC and 
SB-BC) and the symmetry-adapted (SA) one, for the Gaussian and Slater basis sets. 
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symmetry-broken solutions and the symmetry-adapted 
one for the Gaussian basis set of Ref. |^ and the Slater 
basis set of the present study. With the Gaussian basis 
set, the departure of the SB-BC and SB-AC solutions 
from the SA solution occurs for H46 and H54 rings, re- 
spectively. With the Slater basis set, the onset of sym- 
metry breaking takes place for larger rings, i.e. for H50 
and H62 for SB-BC and SB- AC, respectively. In addi- 
tion, for a fixed ring size, the lowering in energy of the 
symmetry-brokeir solutions is smaller with the Slater ba- 
sis set. This is an indication that the Slater basis is better 
than the Gaussian basis, since the amount of symmetry 
breaking is usually larger for poorer wave functions. The 
HF total energies are indeed lower with the Slater ba- 
sis, for example for H42, the (SA) energy per two-atom 
cell is —0.950252 hartree with the Gaussian basis and 
—0.997003 hartree with the Slater basis. As an addi- 
tional verification of the usefulness of Slater functions, 
we can look at Kato's cusp condition [23] at the nuclear 
positions: 



1 |Vp(n.)|c 

2 pirn) 



= Z, 



(3) 



where p(r„) and |Vp(r„)|at). are the density and the 
spherical average of density gradient at the nuclear po- 
sitions, and Z is the nuclear charge. For example, for 
the Hge ring, we find 1.009 for SA, 1.012 for SB-BC, and 
1.034 and 0.980 for SB-AC, close to the ideal value of 
Z = 1. Possibly, symmetry breaking can be further re- 
duced by using a larger Slater basis. 



A. Brief overview of VMC and DMC 

We consider Jastrow- Slater trial wave function of the 
form 



*t(R) = J(R) $(R), 



(4) 



where R designates the electron coordinates, ^(R) is a 
HF determinant and J(R) = e^^^^^^'^'^^'^ is a Jastrow 
correlation factor depending explicitly on the electron- 
electron distances r^ and the nucleus-electron distances 
rji. In VMC, one calculates the energy as the expectation 
value of the Hamiltonian H over the wave function "^t 
by stochastic sampling 



^VMC 



*T(R)i?^T(R)rfR= /"*T(R)^^L(R)rfR 

M 

-^i54R0, (5) 



where ^l(R) = i7*T(R)/*T(R) is the local energy, 
and the M points R; arc sampled from \['t(R)^ by a 
Metropolis algorithm. In DMC, one improves over the 
distribution *I't(R)^ by generating another distribution 
/(R, t) obtained by evolving the importance-sampling 
Schrodinger equation in imaginary time t = it 

-|-/(R,r) = -ivV(R,r) 



-V- /(R,r) 



Wt(R) 
*t(R) 



(i?i(R)-i?T).f(R,r), (6) 



III. QUANTUM MONTE CARLO STUDY 

QMC methods are sometimes believed to produce 
benchmarks in quantum chemistry, approaching the 
electronic correlation problem differently than common 
wave-function-based methods or density-functional the- 
ory. As QMC methods often rely on a HF trial wave 
function, it is interesting to check their sensitivity to HF 
symmetry breaking. We start by giving a brief overview 
of the VMC and DMC methods employed in this work. 



which resembles an ordinary diffusion equation with dif- 
fusion, drift and source terms on the right-hand side. 
This diffusion process is simulated stochastically with 
a population of walkers representing the distribution 
/(R, r). The trial energy Et is adjusted in the course of 
the calculation in order to maintain a stable population 
of walkers. After some iterations, the stationary distribu- 
tion is obtained /(R, r — > 00) = ^fn(R)*t(R) where 
^fn(R) is the fixed- node (FN) wave function, i.e. the 
best approximation to the ground-state wave function 
having the same nodes than the trial wave function. In 
practice, this fixed-node approximation is automatically 



enforced by using ^fn(R-)^t(R-) as a positive probabil- 
ity density, meaning that ^fn(R) niust necessarily be 
of the same sign as \1/t(R-)- The DMC energy is then 
calculated as the statistical average of the local energy 
of the trial wave function over the mixed distribution 
*fn(R)*t(R). 

The nodes of the wave function are the locations of the 
points R where the wave function vanishes. For a sys- 
tem of N electrons in 3 dimensions, they form (3A^ — 1)- 
dimensional hypersurfaccs. A subset of these nodes is 
given by the antisymmetry property of the fermionic 
wave function with respect to the exchange of two elec- 
trons, which implies that the wave function vanishes 
when two same-spin electrons are at the same point in 
space. However, these "Pauli" (or exchange) nodes form 
only {5N — 3)-dimensional hypersurfaccs, and are there- 
fore far from sufficient to determine the full nodal hy- 
persurfaccs (see, e.g., Ref. l23l - [25l for examples of simple 
atomic systems). Likewise, spatial symmetry is gener- 
ally far from sufficient to specify the nodes. For a given 
system, different HF wave functions (of different spatial 
symmetries) share the same Pauli nodes, but otherwise 
generally have very different nodal hypersurfaccs, and 
thus lead to different fixed-node errors in DMC energies. 



B. Computational details 

The QMC calculations have been performed with the 
program CHAMP [23 on a massively parallel IBM Blue- 
Gene architecture using up to 4096 processors. The 
trial wave functions are constructed by multiplying 
the previously obtained HF wave functions by a Jas- 
trow factor consisting of the exponential of the sum of 
electron- nucleus, electron-electron and possibly electron- 
electron- nucleus terms, written as systematic polyno- 
mial and Pade expansions [231 (see also Refs. [2^ [29i ). 
Some Jastrow parameters are fixed by imposing the 
electron-electron cusp condition, and the others are op- 
timized with the linear energy minimization method in 
VMC fljj llSl l30 | , using an accelerated Metropolis algo- 
rithm J3ll 1321 . The orbital and basis exponent param- 
eters are kept fixed in this work. Once the trial wave 
functions have been optimized, we perform DMC calcu- 
lations within the short-time and fixed-node approxima- 
tions (see, e.g., Refs. 1331 - 1371 ). We use an imaginary time 
step of At = 0.01 hartree"^ in an efficient DMC algo- 
rithm having very small time-step errors [3g]. We use a 
target population of 100 walkers per processor, and esti- 
mate statistical uncertainties with blocks of 1000 steps, 
which is larger than the energy autocorrelation time of 
about 50 steps. The statistical uncertainty of the energy 
per H2 cell is smaller than 2 x 10~^ hartree. 

The computational cost of optimizing the two-body 
and three-body terms in the Jastrow factor scales as the 
third power of the number of hydrogen atoms. When 
restricting the Jastrow factor to the two-body terms 
only, the computational cost scales quadratically with 
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FIG. 2: VMC total energies of H4n+2 rings, per H2 cell 
in mhartree, for the symmetry-adapted (SA) and the two 
symmetry-broken (SB- AC and SB-BC) HF solutions from 46 
to 102 hydrogen atoms. The statistical uncertainty is about 
the size of the point. 
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FIG. 3: DMC total energies of H4„+2 rings, per H2 cell 
in mhartree, for the symmetry-adapted (SA) and the two 
symmetry-broken (SB- AC and SB-BC) HF solutions from 46 
to 102 hydrogen atoms. The statistical uncertainty is about 
the size of the point. 



the number of hydrogen atoms, indicating that it is the 
evaluation of the Jastrow factor that dominates the com- 
putational cost and not the evaluation of the Slater de- 
terminant. The large reduction of computational cost 
achieved by removing the three-body terms comes with- 
out a big loss in the VMC energy, and in principle no 
loss at all in the DMC energy. For example, for the H26 
ring system, we find a VMC energy of -13.8894 ±0.0005 
hartree with the three-body term, and — 13.8430±0.0005 
hartree without the three-body term. The computational 
effort is about 20 times more time consuming in the for- 
mer case. We thus use a two-body Jastrow only. Another 
alternative, not employed in this work, is to use a short- 
range three-body Jastrow. 



TABLE II: Total energy and energy differences, per H2 cell in 
mhartree, of the symmetry-adapted (SA) and the symmetry- 
broken (SB-AC and SB-BC) HP solutions, for the Hse ring, 
with the Slater basis set. 
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As the variance of the local energy of a ring of n H2 
molecules is approximately n times the variance for one 
H2 molecule, the statistical uncertainty of the energy 
grows with the square root of n 



ct(E(H2„)) 



ViEL{ii2n)) /ny(Si(H2)) 



M 



M 



(7) 



where M is the number of Monte Carlo iterations. There- 
fore, the statistical uncertainty of the energy per H2 
cell decreases as l/i/n, and thus calculations aiming at 
a given statistical uncertainty for this quantity requite 
fewer steps for increasing ring sizes. 



C. Results 

Figure [2] shows the VMC energies per H2 cell of the 
hydrogen rings with 46 to 102 atoms for the three HF 
solutions. As in the case of M0ller-Plcsset perturbation 
theory and linearized coupled cluster doubles theory, the 
energy ordering of the three solutions is reversed in com- 
parison to HF, the SA wave function giving now the low- 
est total energy and the SB-BC solution giving the high- 
est energy. For Hgg, we show in Table HIl the QMC total 
energies and energy differences of the symmetry-adapted 
and symmetry-broken solutions. For comparison, we also 
report MP2 energies calculated with the same Slater ba- 
sis set. The VMC total energy per H2 cell lie about 50 
mhartree below the MP2 energies, and the energy split- 
tings between the different solutions are smaller than 
those in MP2, which shows the Jastrow factor does a 
good job of describing electron correlation. 

Figure [3] shows the corresponding DMC results. The 
energy ordering is the same as in VMC and MP2, the SA 
wave function giving the lowest DMC total energy, and 
thus the smallest fixed-node error. As shown in Table HIl 
the energy splittings between the different solutions are 
much smaller in DMC. This indicates that DMC is less 
sensitive to symmetry breaking than other correlation 
methods. It is an interesting feature for cases where sym- 
metry breaking cannot be avoided. Of course, symmetry 



breaking would probably be avoided if the orbitals were 
reoptimized in the presence of the Jastrow factor within 
VMC, but this is computationally expensive for large sys- 
tems. 

IV. CONCLUSION 

When HF trial wave functions are used in QMC cal- 
culations, in case of HF instabilities QMC faces the HF 
symmetry dilemma in choosing between the symmetry- 
adapted solution of higher HF energy and symmetry- 
broken solutions of lower HF energies. In this work, 
we have examined the HF symmetry dilemma in hy- 
drogen rings H4„-|_2 which present HF singlet instabil- 
ities for sufficiently large ring sizes. We have shown 
that using a Slater basis set, instead of a Gaussian ba- 
sis set, delays the onset of HF symmetry breaking until 
larger rings and slightly reduces the energy splittings be- 
tween the symmetry-adapted and symmetry-broken wave 
functions. When using these different HF wave func- 
tions in VMC and DMC, we have found that the en- 
ergy ordering is reversed; the symmetry-adapted wave 
function always giving the lowest energy. This con- 
firms previous post-Hartrec-Fock studies in showing that 
these symmetry-broken solutions are bad starting wave 
functions for correlated calculations. The fact that the 
symmetry-adapted wave function gives the lowest DMC 
energy indicates that this wave function has more ac- 
curate nodes than the symmetry-broken wave functions. 
The present experience thus suggests that spatial sym- 
metry is an important criterion for selecting good trial 
wave functions. For systems that are not very large, the 
symmetry-breaking problem could probably be avoided 
altogether by optimizing the orbitals within the quantum 
Monte Carlo calculation, rather than using fixed HF or- 
bitals. 
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